//============================================================================
// Name        : csiena
// Author      : Josh Lospinoso
// Copyright   : MIT License
// Description : Austere, limited implementation of the R package RSiena in
//				 C++. Used for fast maximum likelihood and method of moments
//				 estimation of stochastic actor oriented models (SAOMs)
//				 without any frills or user experience niceties.
//============================================================================

#include "linkEvaluation.h"

namespace csiena {

LinkEvaluation::LinkEvaluation(int actors) {
	this->actors = actors;
	this->expfij = new double[actors];
	this->fij = new double[actors];
	this->pij = new double[actors];
	this->beta = new double[actors];
	this->changeStats = new double[BETA_LENGTH * actors];

	this->pijDenominator=0.0;
	for (int k=0; k<BETA_LENGTH; k++) {
		beta[k] = 0;
	}
	for (int j=0; j<actors;j++) {
		this->fij[j] = 0.0;
		this->expfij[j] = 0.0;
		this->pij[j]=0;
	}
	this->resetChangeStats();
}

LinkEvaluation::LinkEvaluation(const LinkEvaluation &donorEvaluation) {
	this->actors = donorEvaluation.actors;
	this->expfij = new double[actors];
	this->fij = new double[actors];
	this->pij = new double[actors];
	this->beta = new double[actors];
	this->changeStats = new double[BETA_LENGTH * actors];

	this->pijDenominator=0.0;
	// Only difference from regular constructor:
	for (int k=0; k<BETA_LENGTH; k++) {
		this->beta[k] = donorEvaluation.beta[k];
	}
	for (int j=0; j<actors;j++) {
		this->fij[j] = 0.0;
		this->expfij[j] = 0.0;
		this->pij[j]=0;
	}
	this->resetChangeStats();
}

LinkEvaluation::~LinkEvaluation() {
	delete(this->fij);
	delete(this->pij);
	delete(this->beta);
	delete(this->expfij);
	delete(this->changeStats);
}

bool LinkEvaluation::isLinkEvaluation() {
	return true;
}

void LinkEvaluation::preprocess(int i, Digraph* x) {
	if (i < this->actors) {
		this->i = i;
		this->pijDenominator = 0.0;
		this->resetChangeStats();
		for (int j = 0; j<this->actors; j++) {
			this->fij[j] = 0;
			if (i != j) {
				// Outdegree change statistic:
				this->setChangeStats(0, j, 1.0);
				// Reciprocity change statistic:
				if (x->get(j,i)) {
						this->setChangeStats(1, j, 1.0);
				}
				if (BETA_LENGTH > 2) {
					// Log. Transitive triplets change statistic:
					double twoPaths=0.0;
					for (int k=0; k<this->actors; k++) {
						if ( (i != k) & (j !=k)) {
							if (x->get(i,k) & x->get(k,j)) {
								// i -> k, k -> j
								twoPaths++;
							}
							// Uncomment this for generalized two paths
							/*
							else if (x->get(k,i) & x->get(k,j)) {
								// i <- k, k -> j
								twoPaths++;
							} else if (x->get(k,i) & x->get(j,k)) {
								// i <- k, k <- j
								twoPaths++;
							} else if (x->get(i,k) & x->get(k,j)) {
								// i -> k, k <- j
								twoPaths++;
							}
							*/
						}
					}
					this->setChangeStats(2, j, (1+twoPaths) );
				}
				if (BETA_LENGTH > 3) {
					// Log Indegree Popularity Effect
					double alterInPop=0.0;
					for (int k=0; k<this->actors; k++) {
						if ( (i != k) & (j !=k)) {
							if (x->get(k,j)) {
								//  k -> j
								alterInPop++;
							}
						}
					}
					this->setChangeStats(3, j, (1+alterInPop));
				}
				if (BETA_LENGTH > 4) {
					// 3 Cycles change statistic
					double inTwoStars=0.0;
					for (int k=0; k<this->actors; k++) {
						if ( (i != k) & (j !=k)) {
							if (x->get(j,k) & x->get(k,i)) {
								// j -> k, k -> i
								inTwoStars++;
							}
						}
					}
					this->setChangeStats(4, j, (inTwoStars));
				}
				if (x->get(i,j)) {
					this->makeChangeStatsNegative(j);
				}

				// Calculate fij and expfij:
				this->fij[j] = this->beta[0] * getChangeStats(0,j);
				for (int k=1; k<BETA_LENGTH; k++) {
					this->fij[j] += this->beta[k] * getChangeStats(k,j);
				}
			}
			this->expfij[j] = exp(this->fij[j]);
			this->pijDenominator += this->expfij[j];
		}
		// Calculate pij:
		for (int j = 0; j<this->actors; j++) {
			this->pij[j] = this->expfij[j] / this->pijDenominator;
		}
	} else {
		std::cerr << "LinkEvaluation::preprocess: alter index out of bounds: " << i << std::endl;
	}
}

void LinkEvaluation::resetChangeStats() {
	for (int k = 0; k < BETA_LENGTH; k++) {
		for (int j=0; j<this->actors; j++) {
			this->setChangeStats(k,j,0.0);
		}
	}
}

void LinkEvaluation::makeChangeStatsNegative(int alter) {
	if ((alter >= this->actors) | alter < 0) {
		std::cerr << "Evaluation::makeChangeStatsNegative: Actor out of bounds : " << alter << " must be on " << "[0," << this->actors << "]." << std::endl;
	} else {
		for (int k=0; k<BETA_LENGTH; k++) {
			this->setChangeStats(k, alter, -1 * this->getChangeStats(k,alter) );
		}
	}
}

void LinkEvaluation::print() {
	std::cout << "Decision for ego " << i << ": (alter: fij pij changeStatistics)" << std::endl;
	for (int j=0; j<this->actors; j++) {
			printf("%3u : %+6.3f %6.3f", j, fij[j], pij[j]);
		for (int k=0; k<BETA_LENGTH; k++) {
			printf(" |%5.3f", this->getChangeStats(k,j));
		}
		printf("\n");
	}
}
//Getters and setters:

void LinkEvaluation::setChangeStats(int beta, int alter, double value) {
	if ((alter >= this->actors) | alter < 0) {
		std::cerr << "Evaluation::setChangeStats: Actor out of bounds : " << alter << " must be on " << "[0," << this->actors << "]." << std::endl;
	}
	if ((beta >= BETA_LENGTH) | beta < 0) {
		std::cerr << "Evaluation::setChangeStats: Beta out of bounds : " << beta << " must be on " << "[0," << BETA_LENGTH << "]." << std::endl;
	}
	this->changeStats[alter * BETA_LENGTH + beta] = value;
}

double LinkEvaluation::getChangeStats(int beta, int alter) {
	if ((alter >= this->actors) | alter < 0) {
		std::cerr << "Evaluation::getChangeStats: Actor out of bounds : " << alter << " must be on " << "[0," << this->actors << "]." << std::endl;
	}
	if ((beta >= BETA_LENGTH) | beta < 0) {
		std::cerr << "Evaluation::getChangeStats: Beta out of bounds : " << beta << " must be on " << "[0," << BETA_LENGTH << "]." << std::endl;
	}
	double value = this->changeStats[alter * BETA_LENGTH + beta];
	return value;
}

double LinkEvaluation::getFij(int alter) {
	if ((alter >= this->actors) | alter < 0) {
		std::cerr << "Evaluation::getFij: Actor out of bounds : " << alter << " must be on " << "[0," << this->actors << "]." << std::endl;
	}
	return this->fij[alter];
}

double LinkEvaluation::getPij(int alter) {
	if ((alter >= this->actors) | alter < 0) {
		std::cerr << "Evaluation::getPij: Actor out of bounds : " << alter << " must be on " << "[0," << this->actors << "]." << std::endl;
	}
	return this->pij[alter];
}

double LinkEvaluation::getExpFij(int alter) {
	if ((alter >= this->actors) | alter < 0) {
		std::cerr << "Evaluation::getExpFij: Actor out of bounds : " << alter << " must be on " << "[0," << this->actors << "]." << std::endl;
	}
	return this->expfij[alter];
}

double LinkEvaluation::getBeta(int beta) {
	if ((beta >= BETA_LENGTH) | beta < 0) {
		std::cerr << "Evaluation::getBeta: Beta out of bounds : " << beta << " must be on " << "[0," << BETA_LENGTH << "]." << std::endl;
	}
	return this->beta[beta];
}

void LinkEvaluation::setBeta(int beta, double value) {
	if ((beta >= BETA_LENGTH) | beta < 0) {
		std::cerr << "Evaluation::setBeta: Beta out of bounds : " << beta << " must be on " << "[0," << BETA_LENGTH << "]." << std::endl;
	}
	this->beta[beta] = value;
}

void LinkEvaluation::incrementBeta(int beta, double value) {
	if ((beta >= BETA_LENGTH) | beta < 0) {
		std::cerr << "Evaluation::incrementBeta: Beta out of bounds : " << beta << " must be on " << "[0," << BETA_LENGTH << "]." << std::endl;
	}
	this->beta[beta] += value;
}

double LinkEvaluation::getPijDenominator() {
	return this->pijDenominator;
}

int LinkEvaluation::getEgoIndex() {
	return this->i;
}

int LinkEvaluation::getActors() {
	return this->actors;
}

//CUDA Methods here
#ifdef CUDA
void LinkEvaluation::checkForDevices() {
	cudaDeviceProp properties;
	int count = -1;
	cudaGetDeviceCount(&count);
	std::cout << "Found " << count << " CUDA enabled devices." << std::endl;
}
#endif
}
